Source code for hysop.backend.host.python.operator.solenoidal_projection

# Copyright (c) HySoP 2011-2024
#
# This file is part of HySoP software.
# See "https://particle_methods.gricad-pages.univ-grenoble-alpes.fr/hysop-doc/"
# for further info.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.


import functools

from hysop.tools.htypes import check_instance, first_not_None
from hysop.tools.decorators import debug
from hysop.tools.numpywrappers import npw
from hysop.backend.host.host_operator import HostOperator, OpenClMappable
from hysop.core.graph.graph import op_apply
from hysop.operator.base.solenoidal_projection import SolenoidalProjectionOperatorBase
from hysop import __DEFAULT_NUMBA_TARGET__


[docs] def build_projection_filter(FIN, FOUT, K, KK, target=__DEFAULT_NUMBA_TARGET__): assert len(FIN) == len(FOUT) == 3 assert len(K) == len(KK) == 9 args = FIN + K + KK + FOUT try: import numba as nb from hysop.tools.numba_utils import make_numba_signature, prange signature, _ = make_numba_signature(*args) layout = "(m,n,p),(m,n,p),(m,n,p), " layout += "(m),(n),(p), (m),(n),(p), (m),(n),(p), " layout += "(m),(n),(p), (m),(n),(p), (m),(n),(p) " layout += "-> (m,n,p),(m,n,p),(m,n,p)" @nb.guvectorize([signature], layout, target=target, nopython=True, cache=True) def filter_projection_3d( Fin0, Fin1, Fin2, K00, K01, K02, K10, K11, K12, K20, K21, K22, KK00, KK01, KK02, KK10, KK11, KK12, KK20, KK21, KK22, Fout0, Fout1, Fout2, ): for i in prange(0, Fin0.shape[0]): for j in prange(0, Fin0.shape[1]): for k in range(0, Fin0.shape[2]): F0 = Fin0[i, j, k] F1 = Fin1[i, j, k] F2 = Fin2[i, j, k] G0 = K00[i] * F0 G1 = K11[j] * F1 G2 = K22[k] * F2 L0 = KK00[i] * F0 L1 = KK11[j] * F1 L2 = KK22[k] * F2 if (i != 0) or (j != 0) or (k != 0): C0 = (L0 + K10[i] * G1 + K20[i] * G2) / ( KK00[i] + KK01[j] + KK02[k] ) C1 = (K01[j] * G0 + L1 + K21[j] * G2) / ( KK10[i] + KK11[j] + KK12[k] ) C2 = (K02[k] * G0 + K12[k] * G1 + L2) / ( KK20[i] + KK21[j] + KK22[k] ) else: C0 = 0 C1 = 0 C2 = 0 Fout0[i, j, k] = Fin0[i, j, k] - C0 Fout1[i, j, k] = Fin1[i, j, k] - C1 Fout2[i, j, k] = Fin2[i, j, k] - C2 return functools.partial(filter_projection_3d, *args) except ModuleNotFoundError: def filter_projection_3d( Fin0, Fin1, Fin2, K00, K01, K02, K10, K11, K12, K20, K21, K22, KK00, KK01, KK02, KK10, KK11, KK12, KK20, KK21, KK22, Fout0, Fout1, Fout2, ): for i in range(0, Fin0.shape[0]): for j in range(0, Fin0.shape[1]): for k in range(0, Fin0.shape[2]): F0 = Fin0[i, j, k] F1 = Fin1[i, j, k] F2 = Fin2[i, j, k] G0 = K00[i] * F0 G1 = K11[j] * F1 G2 = K22[k] * F2 L0 = KK00[i] * F0 L1 = KK11[j] * F1 L2 = KK22[k] * F2 if (i != 0) or (j != 0) or (k != 0): C0 = (L0 + K10[i] * G1 + K20[i] * G2) / ( KK00[i] + KK01[j] + KK02[k] ) C1 = (K01[j] * G0 + L1 + K21[j] * G2) / ( KK10[i] + KK11[j] + KK12[k] ) C2 = (K02[k] * G0 + K12[k] * G1 + L2) / ( KK20[i] + KK21[j] + KK22[k] ) else: C0 = 0 C1 = 0 C2 = 0 Fout0[i, j, k] = Fin0[i, j, k] - C0 Fout1[i, j, k] = Fin1[i, j, k] - C1 Fout2[i, j, k] = Fin2[i, j, k] - C2 return functools.partial(filter_projection_3d, *args)
[docs] def make_div_filter(target, *args): try: import numba as nb from hysop.tools.numba_utils import make_numba_signature, prange from hysop import __DEFAULT_NUMBA_TARGET__ signature, _ = make_numba_signature(*args) layout = "(m,n,p), (m,n,p), (m,n,p), (m),(n),(p) -> (m,n,p)" @nb.guvectorize([signature], layout, target=target, nopython=True, cache=True) def compute_div_3d(Fin0, Fin1, Fin2, K0, K1, K2, Fout): for i in prange(0, Fin0.shape[0]): for j in prange(0, Fin0.shape[1]): for k in range(0, Fin0.shape[2]): Fout[i, j, k] = ( K0[i] * Fin0[i, j, k] + K1[j] * Fin1[i, j, k] + K2[k] * Fin2[i, j, k] ) except ModuleNotFoundError: def compute_div_3d(Fin0, Fin1, Fin2, K0, K1, K2, Fout): for i in range(0, Fin0.shape[0]): for j in range(0, Fin0.shape[1]): for k in range(0, Fin0.shape[2]): Fout[i, j, k] = ( K0[i] * Fin0[i, j, k] + K1[j] * Fin1[i, j, k] + K2[k] * Fin2[i, j, k] ) return functools.partial(compute_div_3d, *args)
[docs] class PythonSolenoidalProjection( SolenoidalProjectionOperatorBase, OpenClMappable, HostOperator ): """ Solves solenoidal projection (project a 3d field F such that div(F)=0), using spectral methods. """
[docs] def build_divergence_filters(self, target=__DEFAULT_NUMBA_TARGET__): if self.compute_divFin: FIN, K, DIV_IN = self.FIN, self.K, self.DIV_IN K = (K[0], K[4], K[8]) args = FIN + K + DIV_IN self.pre_filter_div = make_div_filter(target, *args) else: self.pre_filter_div = None if self.compute_divFout: FOUT, K, DIV_OUT = self.FOUT, self.K, self.DIV_OUT K = (K[0], K[4], K[8]) args = FOUT + K + DIV_OUT self.post_filter_div = make_div_filter(target, *args) else: self.post_filter_div = None
[docs] @debug def setup(self, work): super().setup(work=work) FIN, FOUT = self.FIN, self.FOUT K, KK = self.K, self.KK self.filter_projection = build_projection_filter(FIN, FOUT, K, KK) self.build_divergence_filters()
@op_apply def apply(self, simulation=None, **kwds): super().apply(**kwds) with self.map_objects_to_host(): for Ft in self.forward_transforms: Ft() if self.compute_divFin: self.pre_filter_div() self.backward_divFin_transform() self.filter_projection() if self.compute_divFout: self.post_filter_div() self.backward_divFout_transform() for Bt in self.backward_transforms: Bt() if self.compute_divFin: self.ddivFin.exchange_ghosts() if self.compute_divFout: self.ddivFout.exchange_ghosts() for Bt in self.backward_transforms: Bt.dfield.exchange_ghosts()